{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "id": "4yPUsdJxSXFq"
   },
   "outputs": [],
   "source": [
    "# @title Copyright 2021 The Cirq Developers\n",
    "# Licensed under the Apache License, Version 2.0 (the \"License\");\n",
    "# you may not use this file except in compliance with the License.\n",
    "# You may obtain a copy of the License at\n",
    "#\n",
    "# https://www.apache.org/licenses/LICENSE-2.0\n",
    "#\n",
    "# Unless required by applicable law or agreed to in writing, software\n",
    "# distributed under the License is distributed on an \"AS IS\" BASIS,\n",
    "# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
    "# See the License for the specific language governing permissions and\n",
    "# limitations under the License."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zC1qlUJoSXhm"
   },
   "source": [
    "<table class=\"tfo-notebook-buttons\" align=\"left\">\n",
    "  <td>\n",
    "    <a target=\"_blank\" href=\"https://quantumai.google/cirq/noise/qcvv/xeb_theory\"><img src=\"https://quantumai.google/site-assets/images/buttons/quantumai_logo_1x.png\" />View on QuantumAI</a>\n",
    "  </td>\n",
    "  <td>\n",
    "    <a target=\"_blank\" href=\"https://colab.research.google.com/github/quantumlib/Cirq/blob/main/docs/noise/qcvv/xeb_theory.ipynb\"><img src=\"https://quantumai.google/site-assets/images/buttons/colab_logo_1x.png\" />Run in Google Colab</a>\n",
    "  </td>\n",
    "  <td>\n",
    "    <a target=\"_blank\" href=\"https://github.com/quantumlib/Cirq/blob/main/docs/noise/qcvv/xeb_theory.ipynb\"><img src=\"https://quantumai.google/site-assets/images/buttons/github_logo_1x.png\" />View source on GitHub</a>\n",
    "  </td>\n",
    "  <td>\n",
    "    <a href=\"https://storage.googleapis.com/tensorflow_docs/Cirq/docs/noise/qcvv/xeb_theory.ipynb\"><img src=\"https://quantumai.google/site-assets/images/buttons/download_icon_1x.png\" />Download notebook</a>\n",
    "  </td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "bd9529db1c0b"
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    import cirq\n",
    "except ImportError:\n",
    "    print(\"installing cirq...\")\n",
    "    !pip install --quiet cirq\n",
    "    import cirq\n",
    "\n",
    "    print(\"installed cirq.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "07034e5e3982"
   },
   "source": [
    "# Cross-Entropy Benchmarking Theory\n",
    "\n",
    "Cross-Entropy Benchmarking (XEB) uses the properties of random quantum programs to determine the fidelity of a wide variety of circuits. When applied to circuits with many qubits, XEB can characterize the performance of a large device. When applied to deep, two-qubit circuits it can be used to accurately characterize a two-qubit interaction potentially leading to better calibration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "1348715511ca"
   },
   "outputs": [],
   "source": [
    "# Standard imports\n",
    "import numpy as np\n",
    "\n",
    "from cirq.contrib.svg import SVGCircuit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "26129d0ff1c0"
   },
   "source": [
    "## The action of random circuits with noise\n",
    "An XEB experiment collects data from the execution of random circuits\n",
    "subject to noise. The effect of applying a random circuit with unitary $U$ is\n",
    "modeled as $U$ followed by a depolarizing channel. The result is that the\n",
    "initial state $|𝜓⟩$ is mapped to a density matrix $ρ_U$ as follows:\n",
    "\n",
    "$$\n",
    "    |𝜓⟩ → ρ_U = f |𝜓_U⟩⟨𝜓_U| + (1 - f) I / D\n",
    "$$\n",
    "\n",
    "where $|𝜓_U⟩ = U|𝜓⟩$, $D$ is the dimension of the Hilbert space, $I / D$ is the\n",
    "maximally mixed state, and $f$ is the fidelity with which the circuit is\n",
    "applied.\n",
    "\n",
    "For this model to be accurate, we require $U$ to be a random circuit that scrambles errors. In practice, we use a particular circuit ansatz consisting of random single-qubit rotations interleaved with entangling gates."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "d940bfde9209"
   },
   "source": [
    "### Possible single-qubit rotations\n",
    "Geometrically, we choose 8 axes in the XY plane to perform a quarter-turn (pi/2 rotation) around. This is followed by a rotation around the Z axis of 8 different magnitudes.\n",
    "\n",
    "These 8*8 possible rotations are chosen randomly when constructing the circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "bb896019c42a"
   },
   "outputs": [],
   "source": [
    "exponents = np.linspace(0, 7 / 4, 8)\n",
    "exponents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "81e2ce9562a5"
   },
   "outputs": [],
   "source": [
    "import itertools\n",
    "\n",
    "SINGLE_QUBIT_GATES = [\n",
    "    cirq.PhasedXZGate(x_exponent=0.5, z_exponent=z, axis_phase_exponent=a)\n",
    "    for a, z in itertools.product(exponents, repeat=2)\n",
    "]\n",
    "SINGLE_QUBIT_GATES[:10], '...'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "72ff411420ef"
   },
   "source": [
    "### Random circuit\n",
    "\n",
    "We use `cirq.experiments.random_quantum_circuit_generation.random_rotations_between_two_qubit_circuit` to generate a random two-qubit circuit. Note that we provide the possible single-qubit rotations from above and declare that our two-qubit operation is the $\\sqrt{i\\mathrm{SWAP}}$ gate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "50f3e9622ff8"
   },
   "outputs": [],
   "source": [
    "import cirq_google as cg\n",
    "from cirq.experiments import random_quantum_circuit_generation as rqcg\n",
    "\n",
    "q0, q1 = cirq.LineQubit.range(2)\n",
    "circuit = rqcg.random_rotations_between_two_qubit_circuit(\n",
    "    q0,\n",
    "    q1,\n",
    "    depth=4,\n",
    "    two_qubit_op_factory=lambda a, b, _: cirq.SQRT_ISWAP(a, b),\n",
    "    single_qubit_gates=SINGLE_QUBIT_GATES,\n",
    ")\n",
    "SVGCircuit(circuit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "b422486e30c9"
   },
   "source": [
    "## Estimating fidelity\n",
    "\n",
    "Let $O_U$ be an observable that is diagonal in the computational\n",
    "basis. Then the expectation value of $O_U$ on $ρ_U$ is given by\n",
    "\n",
    "$$\n",
    "    Tr(ρ_U O_U) = f ⟨𝜓_U|O_U|𝜓_U⟩ + (1 - f) Tr(O_U / D).\n",
    "$$\n",
    "\n",
    "This equation shows how $f$ can be estimated, since $Tr(ρ_U O_U)$ can be\n",
    "estimated from experimental data, and $⟨𝜓_U|O_U|𝜓_U⟩$ and $Tr(O_U / D)$ can be\n",
    "computed.\n",
    "\n",
    "Let $e_U = ⟨𝜓_U|O_U|𝜓_U⟩$, $u_U = Tr(O_U / D)$, and $m_U$ denote the experimental\n",
    "estimate of $Tr(ρ_U O_U)$. We can write the following linear equation (equivalent to the\n",
    "expression above):\n",
    "\n",
    "$$\n",
    "    m_U = f e_U + (1-f) u_U \\\\\n",
    "    m_U - u_U = f (e_U - u_U)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "1cef06bfac12"
   },
   "outputs": [],
   "source": [
    "# Make long circuits (which we will truncate)\n",
    "MAX_DEPTH = 100\n",
    "N_CIRCUITS = 10\n",
    "circuits = [\n",
    "    rqcg.random_rotations_between_two_qubit_circuit(\n",
    "        q0,\n",
    "        q1,\n",
    "        depth=MAX_DEPTH,\n",
    "        two_qubit_op_factory=lambda a, b, _: cirq.SQRT_ISWAP(a, b),\n",
    "        single_qubit_gates=SINGLE_QUBIT_GATES,\n",
    "    )\n",
    "    for _ in range(N_CIRCUITS)\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "9bd38c9d20c8"
   },
   "outputs": [],
   "source": [
    "# We will truncate to these lengths\n",
    "cycle_depths = np.arange(1, MAX_DEPTH + 1, 9)\n",
    "cycle_depths"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "b573f20ea0d2"
   },
   "source": [
    "### Execute circuits\n",
    "Cross-Entropy Benchmarking (XEB) requires sampled bitstrings from the device being benchmarked *as well as* the true probabilities from a noiseless simulation. We find these quantities for all `(cycle_depth, circuit)` permutations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "de9e2414d46f"
   },
   "outputs": [],
   "source": [
    "pure_sim = cirq.Simulator()\n",
    "\n",
    "# Pauli Error. If there is an error, it is either X, Y, or Z\n",
    "# with probability E_PAULI / 3\n",
    "E_PAULI = 5e-3\n",
    "noisy_sim = cirq.DensityMatrixSimulator(noise=cirq.depolarize(E_PAULI))\n",
    "\n",
    "# These two qubit circuits have 2^2 = 4 probabilities\n",
    "DIM = 4\n",
    "\n",
    "records = []\n",
    "for cycle_depth in cycle_depths:\n",
    "    for circuit_i, circuit in enumerate(circuits):\n",
    "\n",
    "        # Truncate the long circuit to the requested cycle_depth\n",
    "        circuit_depth = cycle_depth * 2 + 1\n",
    "        assert circuit_depth <= len(circuit)\n",
    "        trunc_circuit = circuit[:circuit_depth]\n",
    "\n",
    "        # Pure-state simulation\n",
    "        psi = pure_sim.simulate(trunc_circuit).final_state_vector\n",
    "        pure_probs = np.abs(psi) ** 2\n",
    "\n",
    "        # Noisy execution\n",
    "        meas_circuit = trunc_circuit + cirq.measure(q0, q1)\n",
    "        sampled_inds = noisy_sim.sample(meas_circuit, repetitions=10_000).values[:, 0]\n",
    "        sampled_probs = np.bincount(sampled_inds, minlength=DIM) / len(sampled_inds)\n",
    "\n",
    "        # Save the results\n",
    "        records += [\n",
    "            {\n",
    "                'circuit_i': circuit_i,\n",
    "                'cycle_depth': cycle_depth,\n",
    "                'circuit_depth': circuit_depth,\n",
    "                'pure_probs': pure_probs,\n",
    "                'sampled_probs': sampled_probs,\n",
    "            }\n",
    "        ]\n",
    "        print('.', end='', flush=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9902c82e0ff3"
   },
   "source": [
    "## What's the observable\n",
    "\n",
    "What is $O_U$? Let's define it to be the observable that gives the sum of all probabilities, i.e.\n",
    "\n",
    "$$\n",
    "    O_U |x \\rangle = p(x) |x \\rangle\n",
    "$$\n",
    "\n",
    "for any bitstring $x$. We can use this to derive expressions for our quantities of interest.\n",
    "\n",
    "$$\n",
    "e_U = \\langle \\psi_U | O_U | \\psi_U \\rangle \\\\\n",
    "   = \\sum_x a_x^* \\langle x | O_U | x \\rangle a_x \\\\\n",
    "   = \\sum_x p(x) \\langle x | O_U | x \\rangle \\\\\n",
    "   = \\sum_x p(x) p(x)\n",
    "$$\n",
    "\n",
    "$e_U$ is simply the sum of squared ideal probabilities. $u_U$ is a normalizing factor that only depends on the operator. Since this operator has the true probabilities in the definition, they show up here anyways.\n",
    "\n",
    "$$\n",
    "u_U = \\mathrm{Tr}[O_U / D] \\\\\n",
    "    = 1/D \\sum_x \\langle x | O_U | x \\rangle \\\\\n",
    "    = 1/D \\sum_x p(x)\n",
    "$$\n",
    "\n",
    "For the measured values, we use the definition of an expectation value\n",
    "$$\n",
    "\\langle f(x) \\rangle_\\rho = \\sum_x p(x) f(x)\n",
    "$$\n",
    "It becomes notationally confusing because remember: our operator on basis states returns the ideal probability of that basis state $p(x)$. The probability of observing a measured basis state is estimated from samples and denoted $p_\\mathrm{est}(x)$ here.\n",
    "\n",
    "$$\n",
    "m_U = \\mathrm{Tr}[\\rho_U O_U] \\\\\n",
    "    = \\langle O_U \\rangle_{\\rho_U} = \\sum_{x} p_\\mathrm{est}(x) p(x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "9770fc8cf5ba"
   },
   "outputs": [],
   "source": [
    "for record in records:\n",
    "    e_u = np.sum(record['pure_probs'] ** 2)\n",
    "    u_u = np.sum(record['pure_probs']) / DIM\n",
    "    m_u = np.sum(record['pure_probs'] * record['sampled_probs'])\n",
    "    record.update(e_u=e_u, u_u=u_u, m_u=m_u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "e139a1abca2b"
   },
   "source": [
    "Remember:\n",
    "\n",
    "$$\n",
    "    m_U - u_U = f (e_U - u_U)\n",
    "$$\n",
    "\n",
    "We estimate f by performing least squares\n",
    "minimization of the sum of squared residuals\n",
    "\n",
    "$$\n",
    "    \\sum_U \\left(f (e_U - u_U) - (m_U - u_U)\\right)^2\n",
    "$$\n",
    "\n",
    "over different random circuits. The solution to the\n",
    "least squares problem is given by\n",
    "\n",
    "$$\n",
    "    f = (∑_U (m_U - u_U) * (e_U - u_U)) / (∑_U (e_U - u_U)^2)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "2698b1ce5218"
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.DataFrame(records)\n",
    "df['y'] = df['m_u'] - df['u_u']\n",
    "df['x'] = df['e_u'] - df['u_u']\n",
    "\n",
    "df['numerator'] = df['x'] * df['y']\n",
    "df['denominator'] = df['x'] ** 2\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "f526271c8364"
   },
   "source": [
    "### Fit\n",
    "\n",
    "We'll plot the linear relationship and least-squares fit while we transform the raw DataFrame into one containing fidelities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "705fe27d592f"
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# Color by cycle depth\n",
    "import seaborn as sns\n",
    "\n",
    "colors = sns.cubehelix_palette(n_colors=len(cycle_depths))\n",
    "colors = {k: colors[i] for i, k in enumerate(cycle_depths)}\n",
    "\n",
    "_lines = []\n",
    "\n",
    "\n",
    "def per_cycle_depth(df):\n",
    "    fid_lsq = df['numerator'].sum() / df['denominator'].sum()\n",
    "\n",
    "    cycle_depth = df.name\n",
    "    xx = np.linspace(0, df['x'].max())\n",
    "    (l,) = plt.plot(xx, fid_lsq * xx, color=colors[cycle_depth])\n",
    "    plt.scatter(df['x'], df['y'], color=colors[cycle_depth])\n",
    "\n",
    "    global _lines\n",
    "    _lines += [l]  # for legend\n",
    "    return pd.Series({'fidelity': fid_lsq})\n",
    "\n",
    "\n",
    "fids = df.groupby('cycle_depth').apply(per_cycle_depth).reset_index()\n",
    "plt.xlabel(r'$e_U - u_U$', fontsize=18)\n",
    "plt.ylabel(r'$m_U - u_U$', fontsize=18)\n",
    "_lines = np.asarray(_lines)\n",
    "plt.legend(_lines[[0, -1]], cycle_depths[[0, -1]], loc='best', title='Cycle depth')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9703fbf361fd"
   },
   "source": [
    "### Fidelities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "dcb216997aeb"
   },
   "outputs": [],
   "source": [
    "plt.plot(fids['cycle_depth'], fids['fidelity'], marker='o', label='Least Squares')\n",
    "\n",
    "xx = np.linspace(0, fids['cycle_depth'].max())\n",
    "\n",
    "# In XEB, we extract the depolarizing fidelity, which is\n",
    "# related to (but not equal to) the Pauli error.\n",
    "# For the latter, an error involves doing X, Y, or Z with E_PAULI/3\n",
    "# but for the former, an error involves doing I, X, Y, or Z with e_depol/4\n",
    "e_depol = E_PAULI / (1 - 1 / DIM**2)\n",
    "\n",
    "# The additional factor of four in the exponent is because each layer\n",
    "# involves two moments of two qubits (so each layer has four applications\n",
    "# of a single-qubit single-moment depolarizing channel).\n",
    "plt.plot(xx, (1 - e_depol) ** (4 * xx), label=r'$(1-\\mathrm{e\\_depol})^{4d}$')\n",
    "\n",
    "plt.ylabel('Circuit fidelity', fontsize=18)\n",
    "plt.xlabel('Cycle Depth $d$', fontsize=18)\n",
    "plt.legend(loc='best')\n",
    "plt.yscale('log')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "e931726da2af"
   },
   "outputs": [],
   "source": [
    "from cirq.experiments.xeb_fitting import fit_exponential_decays\n",
    "\n",
    "# Ordinarily, we'd use this function to fit curves for multiple pairs.\n",
    "# We add our qubit pair as a column.\n",
    "fids['pair'] = [(q0, q1)] * len(fids)\n",
    "\n",
    "fit_df = fit_exponential_decays(fids)\n",
    "fit_row = fit_df.iloc[0]\n",
    "print(f\"Noise model fidelity: {(1-e_depol)**4:.3e}\")\n",
    "print(f\"XEB layer fidelity:   {fit_row['layer_fid']:.3e} +- {fit_row['layer_fid_std']:.2e}\")"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "xeb_theory.ipynb",
   "toc_visible": true
  },
  "kernelspec": {
   "display_name": "Python 3",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
